source("funkce.R")

skrych = read.table("soubory.dat\\GB-bulk.txt")
sstab = read.table("soubory.dat\\GB-stable.txt")
snestab = read.table("soubory.dat\\GB-unstable.txt")

mkrych = read.table("soubory.dat\\GB-bulk_mat.txt") * 3.65
mstab = read.table("soubory.dat\\GB-stable_mat.txt")
mnestab = read.table("soubory.dat\\GB-unstable_mat.txt")

rownames(skrych) = c("Si1", paste("Ni", 2:4, sep = ""))
rownames(sstab) = c(paste("Si", 1:16, sep = ""), paste("Ni", 17:64, sep = ""))
rownames(snestab) = c(paste("Si", 1:16, sep = ""), paste("Ni", 17:64, sep = ""))


sstab = as.matrix(sstab) %*% as.matrix(mstab)
snestab = as.matrix(snestab) %*% as.matrix(mnestab)

# dvoj1 -------------------------------------------------------------------

vsn = data.frame(a=rep(0,nrow(sstab)), v=rep(Inf, nrow(sstab)))
vns = data.frame(a=rep(0, nrow(snestab)), v=rep(Inf, nrow(snestab)))
for (i in 1:nrow(sstab)){
  for (j in 1:nrow(snestab)){
    k = sqrt(sum((sstab[i,] - snestab[j,])^2))
    if (vsn[i,2] > k) vsn[i,] = c(j, k)
    if (vns[j,2] > k) vns[j,] = c(i, k)
  }
}

dvoj1 = data.frame(NULL)
for (i in 1:nrow(sstab)){
  if (vns[vsn[i,1],1] == i){
    dvoj1 = rbind(dvoj1, c(s = i,vsn[i,]))
  }
}

# dvoj2 -------------------------------------------------------------------

s = unique(sstab[,3])
n = unique(snestab[,3])
for (i in 1 :nrow(sstab)){
  if (sstab[i,3] == 0) sstab[i,3] = s[2] else sstab[i,3] = s[1]
}

vsn = data.frame(a=rep(0,nrow(sstab)), v=rep(Inf, nrow(sstab)))
vns = data.frame(a=rep(0, nrow(snestab)), v=rep(Inf, nrow(snestab)))
for (i in 1:nrow(sstab)){
  for (j in 1:nrow(snestab)){
    k = sqrt(sum((sstab[i,] - snestab[j,])^2))
    if (vsn[i,2] > k) vsn[i,] = c(j, k)
    if (vns[j,2] > k) vns[j,] = c(i, k)
  }
}

dvoj2 = data.frame(NULL)
for (i in 1:nrow(sstab)){
  if (vns[vsn[i,1],1] == i){
    dvoj2 = rbind(dvoj2, c(s = i,vsn[i,]))
  }
}

mean(dvoj1[,3])
mean(dvoj2[,3])

dvoj = dvoj2[,-3]


# anal?za -----------------------------------------------------------------

Dstab = read.csv("ulozeno\\vse_dos_stab.csv",header=T, sep=",", dec=".")[,-1]
Dnestab = read.csv("ulozeno\\vse_dos_nestab.csv",header=T, sep=",", dec=".")[,-1]

Cstab = read.csv("ulozeno\\vse_stab.csv",header=T, sep=",", dec=".")[,-1]
Cnestab = read.csv("ulozeno\\vse_nestab.csv",header=T, sep=";", dec=",")

# snep = (1:nrow(sstab))[!1:nrow(sstab) %in% dvoj[,1]]
# nnep = (1:nrow(snestab))[!1:nrow(snestab) %in% dvoj[,2]]

Ds = Dstab[dvoj[,1], c(-1,-5)]
Dn = Dnestab[dvoj[,2], c(-1, -5)]
r = Ds - Dn

plot(r$kor)
plot(r$d)
plot(r$ddtw)

a = which(dvoj[,1] %in% c(26,49,53))

plot((1:60)[-a], abs(r$kor)[-a])
points(a, abs(r$kor)[a], pch = 4, col = "red")

plot((1:60)[-a], abs(r$d)[-a])
points(a, abs(r$d)[a], pch = 4, col = "red")

plot((1:60)[-a], abs(r$ddtw)[-a])
points(a, abs(r$ddtw)[a], pch = 4, col = "red")
####
as = as.character(Dstab$t[dvoj[,1]])
an = as.character(Dnestab$t[dvoj[,2]])

Cs = Cstab[(Cstab$x %in% as) & (Cstab$y %in% as),]

Cn = data.frame(NULL)
m=c()
for (i in 1:nrow(Cs)){
  x = an[as == Cs$x[i]]
  y = an[as == Cs$y[i]]
  Cn = rbind(Cn, Cnestab[(Cnestab$x == x & Cnestab$y == y) | (Cnestab$x == y & Cnestab$y == x),])
  if(nrow(Cnestab[(Cnestab$x == x & Cnestab$y == y) | (Cnestab$x == y & Cnestab$y == x),]) == 0) m=c(m,i)
}

Cs = Cs[-m, ]

Cs$roz = Cs$ocop/Cs$ovor
Cn$roz = Cn$ocop/Cn$ovor
Cs$roz2 = Cs$ovor/Cs$ocop
Cn$roz2 = Cn$ovor/Cn$ocop

r = data.frame(Cs[,c(1,2)], Cs[,-c(1,2)] - Cn[,-c(1,2)])

bar = c(rep(1,149),2,1,1,2,rep(1,262-153))

plot(abs(r$kor), ylab = "korelace")
r[abs(r$kor)>0.4 &abs(r$kor)<0.5,]
plot(abs(r$d), ylab ="mindist")
r[abs(r$d)>0.035,]
plot(abs(r$ddtw), ylab = "dtw")
r[abs(r$ddtw)>0.015,]
plot(abs(r$out), ylab = "DOHPHS")
r[abs(r$out)>9,]
plot(abs(r$ovor)[r$ovor>-12], ylab = "voronii")
r[abs(r$ovor)>2,]
plot(abs(r$roz[abs(r$roz)<500 | (r$roz %in% c(-Inf,Inf, NA, NaN))]), ylab = "podil")
r[abs(r$roz)>20 & abs(r$roz)<500 & !(r$roz %in% c(-Inf,Inf, NA, NaN)),]
plot(abs(r$roz2))
r[abs(r$roz2)<2 & abs(r$roz2)>1,]


plot((1:262)[-c(150,153)], abs(r$kor)[-c(150,153)], xlab = "Index", ylab = "korelace")
points(c(150,153), abs(r$d)[c(150,153)], pch = 4, col = "red")

plot((1:262)[-c(150,153)], abs(r$d)[-c(150,153)], xlab = "Index", ylab = "mindist")
points(c(150,153), abs(r$d)[c(150,153)], pch = 4, col = "red")

plot((1:262)[-c(150,153)], abs(r$ddtw)[-c(150,153)], xlab = "Index", ylab = "dtw")
points(c(150,153), abs(r$ddtw)[c(150,153)], pch = 4, col = "red")

plot((1:262)[-c(150,153)], abs(r$out)[-c(150,153)], xlab = "Index", ylab = "DOHPHV")
points(c(150,153), abs(r$out)[c(150,153)], pch = 4, col = "red")

plot((1:262)[-c(150,153, which(abs(r$ovor) > 10))], abs(r$ovor)[-c(150,153, which(abs(r$ovor) > 10))], xlab = "Index", ylab = "voronoi")
points(c(150,153), abs(r$ovor)[c(150,153)], pch = 4, col = "red")

plot((1:262)[-c(150,153)], abs(r$roz)[-c(150,153)])
points(c(150,153), abs(r$roz)[c(150,153)], pch = 4, col = "red")

plot((1:262)[-c(150,153)], abs(r$roz2)[-c(150,153)], xlab = "Index", ylab = "podil")
points(c(150,153), abs(r$roz2)[c(150,153)], pch = 4, col = "red")

plot((1:262)[-c(150,153)], abs(r$ocop)[-c(150,153)])
points(c(150,153), abs(r$ocop)[c(150,153)], pch = 4, col = "red")

r2 = r
r2[r2$roz > 40,"roz"] =NaN
r2[r2$roz2 > 4,"roz"] =NaN
rn = data.frame(normalize(abs(r2[,3])))
for (i in 4:10){
  rn = cbind(rn, normalize(abs(r2[,i])))
}
colnames(rn) = colnames(r)[3:10]

celk = rowSums(rn, na.rm = TRUE)

plot((1:262)[-c(150,153)], abs(celk)[-c(150,153)], xlab = "Index", ylab = "celkem")
points(c(150,153), abs(celk)[c(150,153)], pch = 4, col = "red")

